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ABSTRACT 


A few years ago, the class of Essentially Non-Oscillatory Schemes for the numerical sim- 
ulation of hyperbolic equations and systems was constructed. Since then, some extensions 
have been made to multidimensional simulations of compressible flows, mainly in the context 
of very regular structured meshes. In this paper, we first recall and improve the results of an 
earlier paper about non-oscillatory reconstruction on unstructured meshes, emphasising the 
effective calculation of the reconstruction. Then we describe a class of numerical schemes 
on unstructured meshes and give some applications for its third order version. I his demon- 
strates that a higher order of accuracy is indeed obtained, even on very irregular meshes. 


‘Research was supported by the National Aeronautics and Space Administration under NASA Contract 
No. NAS 1-1 9480 while the author was in residence at the Institute for Computer Applications in Science 
and Engineering(ICASE), NASA Langley Research Center, Hampton, VA 23681-0001. 
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1 Introduction 


During the past few years, a growing interest has emerged in building high order accurate 
and robust schemes for compressible flow simulation. One of the difficulties is the appearance 
of possibly strong discontinuities that may interact together, even for smooth initial data. 
One way to avoid this difficulty is to use a Totally Variation Diminishing scheme. Such 
a scheme has the property, at least for ID scalar equations, of not creating new extrema, 
and hence provides a reasonable treatment of discontinuities. TVD schemes have since been 
successfully and widely used with many types of meshes (see for example [1] for a review, 
and among many others [2] for simulations on finite elements type meshes). Nevertheless, 
one of their main weakness is that the order of accuracy falls to first order in regions of 
discontinuity and at extrema, leading to excessive numerical dissipation. 

Various methods have been proposed to overcome this difficulty (for example, mesh 
adaptation [3, 4, 5]), but one promising approach may be the class of Essentially Non- 
Oscillatory schemes (E.N.O., for short), introduced by Harten, Osher and others [6, 7, 8, 
9, 10]. The basic idea of E.N.O. schemes is to use a Lagrange type interpolation with an 
adapted stencil: when a discontinuity is detected, the procedure looks for the region around 
this discontinuity where the function is the smoothest. This reconstruction technique may 
be applied either to the node values [9] or to specially constructed averages over control 
volumes [6, 7, 8]. In the latter case, the approximation is done in a conservative way. This 
enables approximation of any piecewise smooth function to any desired order of accuracy. 

Until now, very few attempts have been made to adapt these ideas to multidimensional 
flows (see for example [9, 1 1]), to smoothly varying grids, or especially to unstructured grids. 
For the latter topic, only preliminary work exists (see [12, 13, 14, 15]). In [14] general 
ideas were presented with a review of the existing ENO methods, but no implementation 
was given. In [15] several algorithms were presented and tested on simple problems (linear 
advection, Burger’s like equations), but their reconstruction algorithms appears to be very 
complicated and costly. There was no study of the numerical stability of the reconstruc- 
tion. Our experience has shown that this point is fundamental. In [12] we studied two 
reconstruction methods based on two different polynomial approximations, and have also 
estimated the behavior of their leading coefficients. This enabled us to design an algorithm 
that has been shown to give third and fourth order approximation. We also discussed the 
choice of candidate stencils: a few stencils suffice. In [13] this reconstruction method was 
implemented for compressible flows problems and tested on a 2D shock tube problem on a 
triangular mesh. In the finite volume scheme used, the control volumes were the triangles 
of the mesh. If one wants this set to be as isotropic as possible, the minimum number of 
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possible stencils is much larger than in the version presented in [12]; its construction is also 
less natural. This makes the scheme of [13] very costly. 

This paper is organized as follows. In section 2 we recall and improve the results ob- 
tained in [12]. These results are valid whatever the type of the underlying mesh (structured 
or not). In particular, we pay much attention to the study of the numerical stability of the 
reconstruction algorithm. In section 3 we present our numerical scheme. In section 4 some 
numerical results are presented. Finally, several general comments are made in the conclu- 
sion. Throughout this paper, the values used are the average values in the given control 
volumes. 

2 The reconstruction problem on unstructured meshes 


Let us first recall basic facts about ID reconstruction. They show why a new method has 
to be introduced for unstructured meshes. We first recall how to interpolate data in an 
essentially non-oscillatory Lagrange fashion, and then how this is used to reconstruct ID 
data. 

Essentially non-oscillatory interpolation. This relies on two well known properties of 
divided differences. Let {y 0 < y\ < . . . < yk) be a stencil, and let [y 0 ? • • * Vk] v be the k + \ th 
divided difference of v y a piecewise smooth real valued function. 

(i) If vis p> n times continuously differentiable in the interval [yo,yk]> then 

[j/ 0 i . ..yk,]v = — + 0{\y k - t/o| n k ) for all k < n and some £ € [yo,yk]- 

(ii) If p < n, admits a jump [u^] in [?/o, Vk], then 

[Vo, ...yk]v = o( ^ L ■) for all k, p + I < k < n. 

\\yh-yo\ v ) 

These relations show that the divided differences remain bounded whatever the mesh size, 
for smooth functions, but go to infinity rather quickly for nonsmooth functions. With the 
help of these two remarks, Harten et. al. have derived the following E.N.O. interpolation 
algorithm: Assume a grid < Vj+i- For any j , first consider S ^ = {yj}- 

CO if |[y>, w+iH < then ^ (1) = {w.yj+iji e]se <?(1) = {w-i.y;}- 
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(ii) Assume that = {y Ja , ■■■, y }k }, a stencil for k+ l t/l order reconstruction, is given (y JO 
ar) d Vj k are the extreme points of the stencil). If J/j 0 , . . . , yk]v\ < \[y J0 , . . . y jk ,y Jk+ i]v 

then S (fc+1) = U{y JO -i}> e ^ se U{j/j*+i} • 

Once the required number of points has been chosen, one can compute Lagrange interpolation 
based on the last stencil of the algorithm: this is the E.N.O. interpolation of v. 

The ID conservative reconstruction. We consider a mesh on IR, {x,}, 6 Z) that may or 
may not be regular. Around each point, x,, we define a control volume [x;_i/ 2 ,x,- + 1 / 2 ] where 
as usual, 

X,' ”1” X l+ 1 

X i+ 1/2 = 2 ' 

Let us consider u, a piecewise smooth real valued function. We let u, denote the average of 
u in [xi_]/ 2 , x )+1 / 2 ]: 

1 f x <+ 1/2 

u, = / u(t ) dt (1) 

x i+1/2 ~ x i-\/2 Jxi-l/2 

There are two classical ways of reconstructing u from its averages: 

(i) Reconstruction via primitive functions. One considers W a primitive of u, say 

w(x) = r u (t) dt. 

Jx\(2 

The values of W at the points Xj +1 / 2 are easily recovered from the data: 

3 

W{x iJr i /2 ) = 5> i+1/a - x,_ 1/2 )u t - 
i=o 

One computes an essentially non-oscillatory reconstruction R(W y n + 1) of W , up to 
the order n + 1, as explained above. Here, one sets y 3 = xj+i/ 2 - The reconstruction of 
if, R(u,n) is defined as: 

, dR(W,n + 1) 

R{u ' n) = Tx 

Clearly the average values of R(u,n) over any [x t *_i/ 2 , ^i+ 1 / 2 ] I s up 

(ii) Reconstruction via deconvolution. Here, the mesh must be regular, x t+ i/ 2 — 
x z _ 1/2 = Ax . One may see equation (l) as the convolution product v of u and the 
characteristic function of [—Ax/2, Ax/2]. One has v(xi) = u,. Then v is reconstructed 
in an essentially non-oscillatory fashion as explained above, where yj = x y Finally, 
one applies a deconvolution operator to 7?(u,n), as in [6] for example, to get R(u,n). 
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Atkins Sz Casper [11] have used a tensor product of lD-reconstructions to derive their 
numerical scheme. This is possible because they assume a regular transformation between a 
Cartesian mesh of [0, 1] x [0, 1] and their computational grid. The same trick has been used 
in Shu et. al. [9]. In the context of unstructured grid, these tricks cannot work for at least 
two reasons: 

(i) The reconstruction via deconvolution needs very regular meshes (each control volume 
can be obtained from any other by translation), 

(ii) The reconstruction via primitive function methods needs to know point values of a 
primitive over any rectangle on the data. Hence, these rectangles must be exactly 
covered by control volumes. As can be seen in Figure 1, this is in general impossible. 

These two remarks show that a straightforward extension of one dimensional ideas is not 
sufficient. 

2.1 Preliminaries 

In the sequel, the symbol IR n [A r , Y] denotes the set of polynomials P in the variables X and 
Y of total degree less than or equal to n : 

P(X,Y) = ± £ 

/=! i+j=l 

The set Y'] is a vector space of dimension N(n) — (»-n)( n + 2 ) ? a basis of which is the 

set of monomials 

{(* - - ^ },+,<* 

where (x 0 ,i/o) is any point of IR 2 . The total degree of P does not depend on the choice 
of (a;o,?/o)- As we will show later, this kind of basis is not the best suited for practical 
calculations. 

Let M. be a mesh of the finite element type. Associated with this mesh, we have a 
triangulation T. We may consider several kind of control volumes, for example the triangles 
of T themselves or the dual mesh (see Figure 2). The dual mesh are constructed as follows: 
for each vertex M , , the control volume is obtained by connecting the midpoints of the edges 
incident on it to the barycenters of the surrounding triangles to which it belongs. Let us 
denote by { C , } the set of control volumes. We only require the following properties: 

• For any i ^ j, C, f) C, is of empty interior, 
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• Ci is connected 


• There is an algebraic dependency of the Cf s on the points of A4 (This is true for the 
two above examples). 

• The boundary of C{ is a locally regular curve (This is also true for the two above 
examples). 

We consider the following problem (problem V or approximation in the mean for short): 

Let u be a regular enough function ( say in L l ). Given N and n two integers t a 
set of control volumes S — {C^} \<i<N, find an element P G IR 7l [X, Y] such that 
for 1 < l < N, 

, , u dx 

< U >Ci = (r 7 —< P >Ci (2) 

1 area[ Cj ( ) ' 


For this problem to have a unique solution, one must fulfill two conditions: 

• N - = N(n) 

• The following Vandermonde matrix must be non singular: 



V = 

[< X'Y* >c.,] , 

+ j<n 




1 

< l < N 


< V > Cll 

< Y >c tl 

... <X n > Cll 

< X n -'Y > Cll • 

.. <V”>c„ 

< X > ClN 

< r >CV 

••• < X n >c iN 




If A s = det V ^ 0, then we say that this stencil is admissible. In that case, there is a unique 
solution to problem V that will be denoted by P u . 

A similar problem was first considered by Barth et. al. [16] for smooth functions, then by 
Harten et. al. [14], Vankeirsblick et. al. [15], and Abgrall [12]. In the three first references 
[16, 14, 15], the authors consider overdetermined systems for two reasons. First, the problem 
V does not always have a unique solution. Second they claim that the condition number 
of the overdetermined system is better than that of problem V. In [12], the same approach 
as here was adopted. To support that choice, one must notice, as explained in remark 1, 
that (3) is generally not singular. Also, the condition number of the linear system depends 
mainly on the basis used for the polynomial expansion, as shown in section 2.4. For these 
two reasons, we have preferred this approach, which also has the advantage of simplifying 
the coding of the global scheme. 
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Remarks: 

(i) We do not know if the admissibility condition admits a geometric interpretation (except 
for n = 1 ). We do not even know whether there is a systematic way of constructing 
admissible stencils, as is the case for the Lagrange interpolation [17]. Nevertheless, 
one may say that, in general, any stencil is admissible: one may consider the equation 
A 5 = 0 as an algebraic surface in IR 2x * for some integer k 2 . This surface is then of 
empty interior, from a topological point of view, so that if <5 is not admissible, one only 
has to change slightly the elements of S for it to become admissible. Nevertheless, the 
condition number of the linear system may be very bad. We will discuss that point in 
section 2.4. 

(ii) This admissibility condition is independent of the basis chosen for expanding the poly- 
nomial P. 

2.2 Some general results about approximation in the mean 

In this section, we give two results on reconstruction in the mean for a function u, which 
may or may not be smooth. These results generalize well known properties of the Lagrange 
interpolation of ID real valued functions, that have been used as a corner stone by Harten 
and his coauthors to design an essentially non-oscillatory reconstruction. We have to empha- 
size that the reconstruction is the mean in not directly related to Lagrange reconstruction. 
Throughout this section, if is an admissible stencil for degree n, the symbol K{S^) 
denotes the convex hull of the union of the elements of S^ n \ 

2.2.1 The case of a smooth function 

In [12] we have shown the following result. Its proof follows easily from Ciarlet & Raviart’s 
proof on Lagrange and Hermite interpolation. 

Theorem 2.1 Let S be an admissible (for degree n) stencil o/IR 2 , let h and p be respectively 
the diameter of K(S) and the supremum of the diameters of the circles contained in K(S). 
Let u be a function that admits everywhere in K(S) an n + I th derivative D n+l u with 

M n+ i = sup{||D" +1 u(a;)||; x G K(S)} < + 00 . 

If P u is the solution of problem V , then for any integer m, 0 < m < n, 

h n+1 

sup{\\D m u(x) - D m P u {x ) ||; z € I<(S)} < CM n+1 — 

2 because we have assumed an algebraic dependency of the control volumes in terms of the points of M 
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for some constant C = C(m,n,5). Moreover , if S' is obtained from S by an affine transfor- 
mation , that is there exists x 0 G IR 2 and an invertible matrix A such that 

C f k G S' iff there exists Ck G S such that C' k = A Ck + %o )? 

then 

This result basically states that if the stencil S is not too flat, i.e. the ratio hj p is not too 
large, then P u will be a good approximation of u. Let us turn now to the case of nonsmooth 
functions. 

2.3 The case of a nonsmooth function 

We only discuss the case of piecewise smooth functions. This class is large enough for our 
purpose. To do the analysis we have to introduce the following property which prevents 
geometrical degeneration: 

Property 2.2 Let us give e > 0. The admissible stencil S ^ belongs to P c n if and only if: 
for any vector U whose components are zeroes and ones and both values are represented y and 
for P the n ih order polynomial defined by 

< P >C{- — Uj , for all Cij G <S (n \ 

the following property holds : the sum of the absolute value of the n-th order coefficieiits of P 
is greater or equal to e } that is 

det(R 0 ■ ■ - Rn(u- i) - " Rr- Rn(u)) ^ 
det(Ro ■ ■ • Rn(u- i) ' * ' Ri ’ * * f^N(n) 

where we have adopted the lexicographic ordering for the monomials {X l Y J }i+j< n , so that 
R k stands for the k th column of the determinant (3) and Ri — U . 

Then we can prove [12] the following theorem that describes the asymptotic behavior of 
the leading coefficients of the approximation in the mean of a piecewise smooth function: 

Theorem 2.3 Let e be a positive real number and S an admissible stencil for degree n such 
that there exists an affine transformation A as in theorem 2.1 for which G P t l - Let 

(£ 0 , 2 / 0 ) be any poirit of the set K{S) and u be a real valued function defined on a open 
subset fi of IR 2 containing We assume that u is in C p_1 , p < n } on and , except 
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on a locally C curve , admits a continuous and bounded p i * 1 derivative with a jump [Z) p u], 

|[D p u]| > M p > 0. Then, the highest degree coefficients of the Taylor expansion of P u satisfy 

E M>C(n,p,c)^ (5) 

i+j=n n 

where C(n,p,e) is a constant independent of S and invariant by affine transformation. 

2.4 Study of the linear problem to solve for the reconstruction 

In this section, we study the numerical system one solves to get P u from the data. We 
consider two kinds of expansions of P u : 

(i) The “natural” expansion: for any point (ar 0 ,jfo) £ IR 2 , 

Pu = ^2 a ij(X — xoy(Y — y 0 y (6) 

i+j<n 

(ii) An expansion using “barycentric” coordinates that we now describe: let } C 2 , C 3 ,...,C h 

be an admissible stencil. Hence, at least one subset of three elements of S ^ is an ad- 
missible stencil for n = 1. We may assume that the set {C u C 2 ,C 3 j is admissible. We 
consider the three polynomials A, of degree 1 defined by: 

< A, > c = S 3 { , 1 < * < 3, 1 < j < 3. (7) 

The symbol Si is the Kronecker symbol. Clearly, we have 

Aj + A 2 T A3 = 1 . 

These polynomials are the barycentric coordinates of the triangle constructed on the 
gravity centers of C\, C 2 , and C 3 . In order to get expansion 6, a strategy may be to 
look first for the expansion of the polynomial P u in terms of power of A 2 and A 3 : 

P = E «oA'A J 3 (8) 

i+j<n 

and then to get the Taylor expansion of P u around the barycenter of C\ from (8) (the 
theorems 2.1 and 2.3 give the behavior of the leading coefficients of P u whatever the 
point chosen in the convex hull of S ). 

In order to get the expansions (6) or (8), one has to solve linear N(n) x N(n) systems: 

B{ci 00 ■ ■ ■ «o„) T = (< U >c,, • • • < u >c, N(n) ) T (9) 

where the matrix B is obtained by taking the average of (X - xq)'{Y - y 0 ) J for (6) and A* 2 A^ 
for (8). Let us now study the properties of these linear systems. 
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2.4.1 Case of expansion (6) 

A very easy consequence of the inequality (5) is that: 

Proposition 2.4 Let us assume that the conditions of theorem 2.3 holds , and let h be the 
supremum of the diameters of the spheres containing /\(«S^). Then the. condition number 
of system (9) is at least 0{h~ n ) for h small enough. 

Proof: For the sake of simplicity we consider the following norm on IR n [A", V']: for P = 
£i+i<„ a tJ (X - x o y(Y - y 0 ) j , ||P|| = Ei+i<n Kl- 0n IRN(,l) ’ we cons i^ er the L l norm. Let 
U be a set of data for the right hand side of (9), and consider the perturbation 8U, 

8U = (0---e---) T 

where e is at the l th position, l > N(n - 1) + 1. All the other entries of U are zero. If one 
considers the function u defined on |JC; by 


x € (~>i , u(x') — 8Ui, 

one can apply theorem 2.3. Hence, the perturbation 8P has a norm satisfying 

ll«p|| > E M > cf 

i+j=n 

since \ \SU\ \ — e. This complete the proof. □ 

This fact is well known for ID Lagrange interpolation and has motivated the search 
for more efficient algorithms, such as the Newton algorithm. There exist algorithms that 
generalizes it [18, 19]. They involve numerous solutions of linear systems, so that we have 
preferred a more classical approach (see section 2.5), for which the coefficients of the linear 
systems are obtained from the a bary centric” coordinate expansion (8) as it is explained now. 


2.4.2 Case of expansion (8) 

In the case of expansion (8), we have the following result: 

Proposition 2.5 If property 2.2 holds for some e > 0, then the condition number of the 
system (9) for the expansion (8) is bounded above and below by constants independent of h } 
the supremum of the diameters of the circles containing K{S ^). 

Proof: The proof is also based on that of theorem 2.3. As in proposition 2.4, the only thing 
that we have to do is to study the effect on the aif s of a perturbation 5U . We denote by P 
the polynomial whose averages are defined by 8U . The proof can be achieved in two stages: 
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(i) Let B any invertible matrix. Consider the stencil 

«SW = {B\Ci, ] +x 0 \ 

l L tjj . u Ji<j</V( n ) 

for any x 0 . It is clear from the definition of the A*’s that A,(x) = A*(x) if x = 
Z?x + Xo- Hence, the sum 5(P) of the absolute values of the coefficients of P in the 
basis A^(x) A^(x) is the same as that of the development of P in the basis A^xJA^x). 
This is a homogeneity property. 

(ii) Since the set of stencils defined by property (2.2) is compact, S(P) is bounded below 
and above, independently of B, hence independently of h: 


Ci >F>C 2 {e } n) >0. 


The constant Ca(e, n) is larger than zero because SU ^ 0 


This achieves the proof □ 


2.5 The explicit calculation of the reconstruction 

From the previous results, the evaluation of the coefficients a tJ in (6) is done through those 
of (8) and hierarchally. For the sake of simplicity, we assume that for any p < n, the set S ^ 
of the N(p) first elements of is an admissible set for order p. This can be achieved with 
a suitable numbering of the elements of The idea is that instead of looking directly 

for the coefficients of P^ n \ to get first those of all of the P^’s, the reconstruction over $<*>, 
for 1 < k < n. Then to construct those of P^ n b In the ENO algorithm described in section 
2.6, this involves no extra cost and simplifies the evaluation of the a^-’s. This also has the 
advantage of reducing the size of the linear systems and also of improving their condition 
numbers. 

Assume that P^ x \ . . . , P ^ are known. We first compute the coefficients of p( p+1 ) — p( p ), 

p(r+i)_p( P ) = Y: a' t] A’A J 3 

*+i<p+ 1 

by solving the linear system 

V. (-) = (/’’ n " ) ( 7 

\ £2 / \ ^ P+l p+l J \ 02 

In equation (10), a\_ (respectively £ 2 ) stands for the coefficients {a,j};+j< p (resp. (aij}t+i=p+i)- 
The block matrices A v , B p p+1 , C p p+ i and D p p+] are defined according to this decomposi- 
tion. In particular, we notice from the hypothesis that A p is invertible. 


Ml 
u 2 


( 10 ) 
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From the conservation property, we get = 0, so that the system (10) can be split: 

Qi = Bp 

O') 

Cp p+1 Ap Bp p+l + Dp p-j-iJ ^2 ^2 

Since <S (p+1) is admissible, E v = [— C p p +iA~'B p P+1 + D p p+1 ] is also invertible, so that one 
can get 02 , then aj_, and finally the coefficients of p( p+1 h 
Simple manipulations show that 

4- 1 _ f p+i^p -1 ^? P+i^p 1 + ^p 1 “^p^p p+i^p 1 ^ 

p+1 ~v -^p-^pp+iV V J 

so that one can quite easily go to the next step. In our case, since the total degree of the 
reconstruction is less or equal to 4, at most two stages of this method are needed. 

Finally, one must notice that the condition number of this method is always better that 
that of the original one because it depends only on part of the original system. 

2*6 The E.N.O* reconstruction 

In [12], we have found that only a few stencils were indeed necessary to achieve an essentially 
non-oscillatory reconstruction of a piecewise smooth function. This set has to be as isotropic 
as possible. Moreover, the ENO reconstruction was found to achieve the expected order of 
accuracy for smooth functions, even on very irregular meshes. In what follows a ZJ always 
stands for any of the coefficients of the reconstruction P in the natural basis, {(X — Xo) 2 (E — 
j/o) J }- Let us describe our procedure up to fourth order: 

(i) Let us start from a given cell, Co assigned to a point of Xf, say (x 0 ,yo)* 

(ii) Consider all the triangles having (#o, Vo) as a vertex, and choose the one, say T m i n , 
that minimize 

J2 K+ 

i+j = I 

Here, <S^ is the set of control volumes corresponding to the vertices of T m i n , (see 
Figure 3-a). For a regular unstructured mesh, there are about six possible triangles. 

(iii) Consider T mm . For any of its three edges, consider the three triangles, Ti, T 2 , T$ as in 
Figure 3-a. There are three possible configurations. We choose the one that minimizes 
the sum 

E M- 

i+j = 2 
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(iv) Consider, as in Figure 3-b, the configuration for third order. It is obtained as fol 
lows: for a stencil made of the control volumes associated with the vertices of 
{T min ,Ti,T2,T 3 }, one may consider its “edges” made of the external sides of {T 2 ,T 3 }, 
{T 2 , T mi n], { Tmin , T 3 }. Consider one of them, say {T 2 , T 3 }, and the vertices «, 0 and 7 . 
Since the triangulation is conformal, there exist a triangle X4 ^ T 2 011 fbe other side of 
[cr, fi\. Similarly, T$ for [/?,7]- Then one can construct Tg, T 7 , Tg and Tg, analogously. 
The stencils for fourth order reconstruction are the union of <S (2) and the control vol- 
umes associated with the additional vertices of either {X4, T5, Tg, TV} or {T4, T5, T9, Tg} 
or {Tj,T 5 ,T 7 ,T 8 }. For a stencil <S (2) , there are at most 12 stencils for fourth order 
reconstruction. 

The situations seem to be become more and more complicated as the degree increases. 
Nevertheless, there is a very easy way to simplify it, so that at each level only three new 
stencils for the n -f \ th order have to be considered over those of a n th order one, just as 
from second order to third order. Assuming a mesh M , we want to derive aH I th order 
reconstruction method. The idea is to work with the control volumes defined for a mesh Ad , 
the points and the triangulation of which are constructed from those of M by adding, for each 
triangle of M, the points and the triangles associated with the P k Lagrange interpolation 
[ 20 ]. 


3 A class of high order numerical scheme for com- 
pressible flow simulations 


3.1 The Euler equations 

Let us quickly recall elementary things about the Euler equation of a calorically perfect gas. 

aw 9 F{W) dG(W) (12) 

dt dx dy 

As usual, in equation (12), W stands for the vector of conserved quantities and F (respec- 
tively G ) is the flux in the x direction (resp. y direction): 


F(W) = 


p u 

p u 2 + p 
p uv 
u(E + p) 


G(W) 


p V 

p uv 

p v 2 + p 

v(E + p) 


with initial and boundary conditions. In equation (13), p is the density, u,v are the com- 
ponents of the velocity, E is the total energy, and p the pressure, related to the conserved 


12 



quantities by the equation of state: 


P = (l - 0 


(14) 


The ratio of specific heats, 7, is kept constant. 

It is well known that the system defined by equations (12), (13) and (14) is hyperbolic: 
for any vector n = (n x ,n y ), the matrix 


„ dF dG 
A!i ~ n ‘dW +n, dW 


(15) 


is diagonalizable and has real eigenvalues and eigenvectors. Let us describe now the con- 
struction of a k th order scheme. 


3.2 Finite volume formulation 

We consider a mesh M. and the control volumes as in figure 2. The semi discrete finite 
volume formulation of (12) is: 

1^(0 = ^ L = Ci(t) (16) 

at area{Li) Jdc , 

Here W (<) is the (spatial) mean value of W (x, t ) at time t over C t , n = (n r , n y ) is the outward 
unit normal to dCi, and Ta = n x F + n y G. We first describe the spatial approximation of 
(16), then the temporal discretization of the resulting set of ordinary differential equations. 
Finally, we detail the boundary conditions. 

3.2.1 Spatial discretization 

For the sake of simplicity we define the integer number p such that either k — 2 p or k = 2 p+ 1 . 
The first step is to discretize C t (t) up to k th order. First, we can rewrite area(Ci)Ci(t) as 

/ Frt[W(x,i)\dl = Y, [ F*[W(x,t)]dl (17) 

JdCi Jr, 

where, as in figure 2, the set IVs is that of the linear edges of C{. On each T s , n is constant. 
We consider on any the p Gaussian points {G , ;}k/< p associated to the Gaussian formula 
of order 2p+ 1. The integral J Ts F,i[W{x,t)]dl is approximated by 

X>e»,,(0 (is) 

;=i 

where term Qa,i{t) is defined now. Set C } is the other control volume of which is a part 
of the boundary. In C, and C } , one computes the ENO reconstructions at time t of W , 
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/?,[VF( . , f),fc] and Rj\W{ . , t),k], up to order Jc. The ENO reconstruction of section 2 is 
applied to the physical variables, then one deduce the conserved ones. From that, we set, in 
equation (18): 

On,,(i) = Jf aman {R i [W( . ,t),k](Gi), Rj[W( . ,t),k](G,)}. (19) 

In equation (19), y r ? ieraann may be any of the available Riemann solvers. In all the example 
below, we have chosen Roe’s Riemann solver with the Harten-Hyman entropy correction. 


3.2.2 Temporal discretization 


The equations (16), (IT), (18) and (19) define a finite set of ordinary differential equations 
that we write as: 

|Wi(0 = £(t) (20) 

In (20), £,-(<) is the discrete version of Cfit). This equation is discretized by the k th order 
version of the Runge-Kutta scheme of Shu [9]: 


' W, {1) = + fii m C ( ”% l = 1 , 2 , 

wj 0) = w?, w} p) = w; i+ \ 


c\ m) = 


( 21 ) 


The order of accuracy, as well as its TVD properties, is achieved by adequate sets of coeffi- 
cients ai m , film, and p (see [9] for details). 


3.2.3 Boundary conditions 


Let T be the boundary of the computational domain and ft be the outward normal unit on 
T. We assume that T is divided into two parts, T = Toff Too, on which different boundary 
conditions will be used. Here, To represents a solid wall, while represents the far-field 
(inflow or outflow). 

We do not treat a boundary conditions by forcing the value of a variable to a prescribed 
boundary value, but consider instead the integral formulation (16) and apply boundary 
condition by modifying the flux integrals on dC t for those cells with Tf)^ / 0. 

For example, for a vertex i located on r 0 , we do not impose the slip condition U • n — 0 
but take this condition into account in the evaluation of the convective flux 



Dae, 


Fn x + Gn y 


0 

Jr 0 p| ac > prix 

fr 0 fl dc, P n y 
0 
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The pressure integrals are computed as: 


/ pn x ~ p, / 

J r 0 fl act Jr 


r 0 f| 9C i 


r x. 


Lf ]x , pn ^ p 'i 


r 0 fl 9(: ' 


n r 


For a vertex located on Too, we again use an approximate Riemann solver. We define a 
far field state W^, n, = / roo p| 9C , " and se U * n agreement with what has been done in the 
interior of the computational domain, 


I Fn x + Gn y = <D(W„ n,). (22) 

T oo f'] dC{ 

In equation (22), is a numerical flux function. For simplicity, we have chosen a modified 
Steger- Warming flux instead of the Roe flux, 


$(W i ,W o0 ,n) = AtW i + AzW c 


The matrices A\ and A ^ are the positive and negative parts of the matrix A a defined in 
(15) and evaluated for W = W,. 

In all the examples we have treated below, the boundaries were either fully subsonic 
or fully supersonic, so that the procedure was really simple, contrary to what would have 
appeared in a mixed type boundary condition. 

Finally, we have reduced the order of accuracy of the reconstruction for cells that are too 
near to the boundary. For them, a proper calculation of the ENO stencil may be impossible 
because the set of possible stencils is biased in one direction due to the boundary. For the 
third order scheme, these cells are those related to a mesh point that belongs to a triangle 
having at least one point on the boundary. For the fourth order scheme, they are those 
belonging to a triangle which shares a vertex with a cell for which a reduction of order must 
be done for third order. 


3.3 Positivity of the density and the pressure 

As pointed out by Harten et. al. [22], in some situation and for a very few cells, the ENO 
reconstruction of the density and pressure may lead to negative values. For these cells, and 
these cells only, following [22], we reduce the order of accuracy with the following inductive 
method (w is either the density or the pressure, to, is its average on C,). Consider, in C,, the 
reconstruction 

£ a r ,(X-x,Y(Y-y,)\ 

f=0 P^rq-l 

If zr = 2 Z P+9= , Ml(x - x g y{y - y g y I > a\w t \ at a Gaussian point ( x,y ), then the recon- 
struction, for that point, is set to R[w,n — l](a:,j/). Then, we repeat the test if necessary. 
Usually, the parameter a is set to 0.95. 
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In all the simulations we have done, the tests were positive for a very small set of cells 
and a zeroth order reconstruction was never used. They were never positive for the second 
order scheme. The number of points causing difficulty is problem dependent. For example, 
only three points caused occasional problems for the facing step problem with a 5000 node 
mesh. These points were all located in the front shock. 

4 Numerical tests 

All the example we propose now have been computed with the second and third ENO 
schemes. The ratio of specific heats, 7 is always set to 1.4. 

4.1 A linear advection problem 

In order to test the precision of these scheme, we have computed the advection of sine wave 
on a sequence of totally unstructured meshes with an increasing number of points. The con- 
vection velocity was parallel to the x axis but since the meshes was totally unstructured, this 
is not a privileged direction. Figure 4 shows, in the abscissa, the logarithm of the maximum 
radius of the circumscribed circles of the triangles of the meshes, and in the ordinate, the 
logarithm of the maximum absolute value of the difference between the computed values. 
The exact solution is also indicated. The slopes —2 and —3 are indicated, so that one can 
see that the expected order of accuracy is indeed achieved. Figure 5 shows, for the medium 
mesh, a cross section of the computed solution for the second and third order scheme. The 
exact cross section is also indicated. One can see that the main differences lies at the extrema 
of the sine, as expected. 

4.2 A Shock tube problem 

We have set up a two dimensional shock tube problem in the square [0, 1 ] x [ 0 , 1 ]. Its 
boundary are solid. The initial conditions are: 

f , = i.o 

forx < 0.5 and | y — 0.5| < 0.25, < u = v = 0.0 

l P = 1.0 

' p = 0.125 

else < u — v — 0.0 

, P- 0.1 

The mesh is completely unstructured with 2127 nodes and 4088 triangles. The velocity field 
obtained by the third order scheme at time t — 0.9 is displayed in Figure 6 . The differences 
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between both results are more clearly visible in the near stagnation zone. In order to better 
represent that area, we have removed from the velocity field all the points for which the sum 
of the absolute values of the two components is larger than 0.15. The result is shown in 
Figures 7 (second order) and 8 (third order). One can clearly observe that the number of 
small structures of the flow is much more important in Fig. 8 than in Fig. 7. The shocks in 
the upper and lower part of the pictures are also resolved differently. Their location is also 
different, though this can be seen only by superimposing the pictures. 

One should also mention that this test is not particularly easy for our method. After a 
short time, the shock reflects from the wall. The reflected shock interacts with the others 
structures of the flow, leading to interactions between the various kind of discontinuities 
and with the smooth parts of the flow. The multiple regions, as shown on our figures, with 
different kind of discontinuities (contact and shock) are resolved by our method without any 
special tricks. 

4.3 A Mach 3 wind tunnel with a step 

We have run this test case, documented in [23], for the second order and third order ENO 
schemes on a 5140 node, 9958 triangle mesh. This discretization corresponds to the medium 
mesh used in [23]. A portion of it is displayed in Figure 9. It is totally unstructured. The 
conditions of the problem are the following: a uniform Mach 3 flow is set in a channel. At 
the initial time, a step of relative height 0.2 is installed in the channel. The channel length 
is 3 and the step is located at 0.6. This situation creates a shock that reflects on the upper 
part of the channel then evolves to a lambda shock as time increases. It interacts with the 
upper part of the step. A weak shock is also created by the expansion wave at the corner. 
This shock interacts with the reflected one creating a slip line. The location of this slip line 
is highly dependent on the boundary conditions set at the corner. 

Here, no special treatment is done, contrary to what was advocated in [23], so that 
the quality of the second reflected shock is poor. We only want to verify the effect of the 
increasing order of accuracy on the solution, so that we will only look at the first reflected 
shock. The solutions of Figure 10 (second order ENO) and Figure 11 (third order ENO) are 
shown. A clear improvement on the thickness of that reflected shock can be seen from the 
horizontal cross section of the density at y = 0.5, Figure 12. The slip line coming from the 
lambda shock is also more visible in Figure 11 than in Figure 10 as well as the weak shock 
near the corner. 


17 


4.4 Reflection of a shock on a wedge 

This problem is also well documented in the literature. In order to achieve a correct solution, 
one has either to use very fine meshes or adapted meshes (see [5] for example). We have 
chosen a case where the planar shock enters from the left in a quiescent fluid. Its Mach 
number is M$ = 5.5 and is defined for the flow values in the quiescent fluid where the 
density is set to 1.4 and the pressure to 1. One expects a double Mach reflection. 

The mesh has only 8569 points and 16806 triangles. A part of it is shown in Figure 13. 
The density contours of the two calculations are displayed in Figures 14 (second order) and 
Figure 15 (third order). A very clear improvement of the slip line coming from the Mach 
stem can be observed. The second triple point can also been observed in Figure 15, though 
it is of poor quality because of the insufficient resolution of the present mesh. It is, however, 
totally obscured in Figure 14. Generally speaking, all the discontinuities are better resolved 
by the third order scheme. 

5 Conclusions 

A third order ENO scheme has been derived for triangular unstructured meshes; this demon- 
strates the possibility of deriving ENO schemes for unstructured meshes. We indicate how 
to build higher order ENO schemes and give some comments on the numerical stability of 
the reconstruction step. 

Our new scheme has been tested on a set of well known test cases and compared to 
a second order scheme. In all cases, the results are clearly improved. Our results also 
demonstrate the scheme’s robustness. The cost of the scheme is four times that of the second 
order scheme (on a Cray YMP) even though the new code is far from being optimized. In 
particular, no optimization has been done in the ENO reconstruction procedure, the most 
expensive routine, so that the factor of four is clearly a poor upper bound on the cost ratio. 

In the near future, we will derive the fourth order version of this class of schemes. The 
two schemes will then be coupled with a dynamic adaptation procedure [3] to improve their 
efficiency. 
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Figure 1: Covering 


of the rectangle [x 0 , xx] x [y 0 , J/i] by triangular control volumes 
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(») (b) 

Figure 3: Stencils for third and fourth order reconstruction 
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Figure 4: Representation of the logarithm of the L error in term of the logarithm of the 
maximum radius of the circumcircles. 
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Figure 6: Shock tube : velocity field at time t = 0.9 
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Figure 9: Portion of the mesh used for the step case 
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Figure 10: Density contours for the second order ENO solution, t=4, min=0.329, max=4.64. 
Density contours from 0.287 to 4.584, A p = 0.14 



Figure 11: Density contours for the third order ENO solution, t=4, min=0.287, max=4.584, 
A p = 0.14 
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Figure 12: Cross-section of the density, y = 0.5 O : second order, + : third order 


32 




Figure 14: Refection of a planar shock by a wedge : density contours, second order solution. 
Min=1.4, Max=17.3. Contour from 1.4 to 19.088, A p= 0.36 
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Min=1.4, Max=19.088, A p = 0.36 
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